function [P_S, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs, educ_mat)
    
    global  a bet0 bet1 pe tau w ub lb N_w N_a sigmae b ...
            thetas w_bar unemp xi eta  var_phi share curve
    
    % Pre-allocate
    u_s = zeros(N_w, N_a, length(Ss), 2);
    F_es = zeros(N_w, N_a, length(Ss));
    wprimeS = zeros(N_w, N_a, length(Ss));
    [~,id]=min(rs);

    for c = 1:length(Ss)
        E_S = (a' .* (bet0 + bet1*Ss(c)).^xi)  .* f_w(w) ;         
        btilde = unemp + b .* f_w(w) + pe .* E_S;
        
        % For use in other programs
        e_S = educ_mat(:,:,c);
        F_es(:,:,c) = eta * (pe + e_S) .* E_S ./  f_w(w);
        wprimeS(:,:,c) = unemp + (b  + F_es(:,:,c)) .* f_w(w)  ;
        
        u_s(:,:, c, 1) = (max((w - rs(c) - tau*e_S.^var_phi).*( wprimeS(:,:,c) ), 0) );
        
        % if wage is too low to afford housing, allocate to cheapest neighborhood
        if (c == length(Ss))
            mysub = sum(u_s, 3:4);
            [mrow, mcol] = ind2sub(size(mysub), find(mysub == 0));
            u_s(mrow,mcol, id,1) = 1e-20;
        end
    end

    % Amenities shock
    for c = 1:length(Ss)
        u_s(:,:,c,2) = u_s(:,:,c,1) .* exp(log(thetas(c)));
    end

    % Calculate choice probabilities
    my_log_diff = (log(u_s) - max(log(u_s), [], 3)) ; 
    numer = exp(my_log_diff .* sigmae ); 
    P_S = numer ./sum(numer , 3);
      
end